Phonons for the FCC Lattice

The Lattice vectors

The direct approach:

Just time them in by hand:

lvecs={{-(a/2),-(a/2),0},{-(a/2),0,-(a/2)},{-(a/2),0,(a/2)},{-(a/2),(a/2),0},{0,-(a/2),-(a/2)},{0,-(a/2),(a/2)},{0,(a/2),-(a/2)},{0,(a/2),(a/2)},{(a/2),-(a/2),0},{(a/2),0,-(a/2)},{(a/2),0,(a/2)},{(a/2),(a/2),0}};
Grid[lvecs,Alignment->"a"]

The Clever Approach

Make a list of all combinations, many of which we don't want:

bigTable=Table[(a/2){i,j,k},{i,-1,1},{j,-1,1},{k,-1,1}]

The above is a list of lists of  lists.   We have to strip it to a list of vectors (a list of lists).

bigList=Flatten[bigTable,2]

Now select only those elements that have the right Cartesian length:

lvecs=Select[bigList, #.# == ((a^2)/2) &]

The Calculation

The crystal momentum:

qvec={qx,qy,qz}

Now just make the matrix :

dynmat=Sum[ lvec=lvecs[[i]];
((Sin[((qvec.lvec)/2)])^2)(A* ((Outer[Times,lvec,lvec])/(lvec.lvec))+B*IdentityMatrix[3] ) ,{i,1,Length[lvecs]}]//Simplify

MatrixForm[dynmat]

We can extract the eigenvalues:

evals=Eigenvalues[dynmat];

Let's plot along the (100) direction.  That means that qy=qz=0.  We will further assume that A=2 and B=1.

evalsPlot=evals/. {A-> 0, qy-> 0, qz-> 0, a-> 1, B-> 1}

This is not very informative!   Let's force Mathematica to write these equations with radicals:

evalsPlot=ToRadicals[evalsPlot]

We see that they are all degenerate.

Plot[Evaluate[(Sqrt[evalsPlot])],{qx, -π, π}]

Or, alternatively, we could set qy=qz=0 in the dynamical matrix:

dtemp=dynmat/. { qy-> 0, qz-> 0} // FullSimplify;
MatrixForm[dtemp]

Now we see how the eigenvalues of the matrix depend upon A and B.   Let's not set A =0 in this case.

temp1=Eigenvalues[dtemp];
temp2=temp1/.{ A-> 1, B-> 1,a-> 1};
Plot[(Sqrt[temp2]),{qx, -π, π}]

Let's look in the 111 direction.

dtemp=dynmat/. { qx-> q,qy-> q, qz-> q} // FullSimplify

temp1=Eigenvalues[dtemp];
temp2=temp1/.{ A-> 1, B-> 1,a-> 1};
Plot[(Sqrt[temp2]),{q, -π, π}]

How about the polarizations?

evecs=Eigenvectors[dtemp]/.{ A-> 1, B-> 1,a-> 1}

ϵ1=Simplify[evecs[[1]]/(Sqrt[evecs[[1]].evecs[[1]]])]
ϵ2=Simplify[evecs[[2]]/(Sqrt[evecs[[2]].evecs[[2]]])]